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Abstract 

We study the propagation of stationary waves in disordered non-linear media 
described by the nonlinear Schrodinger equation and show that for given 
boundary conditions and a given coherent wave incident on the sample, the 
number of solutions of the equation increases exponentially with sample size. 
We also discuss the ballistic case, the sensitivity of the solutions to the change 
of external parameters, the similarity of this problem to the problem of spin 
glasses and time-dependent solutions. 

I. INTRODUCTION 

In this article we consider the propagation of a stationary coherent wave described by a 
field 0(r), in a nonlinear elastically scattering medium. Though we believe that our results 
are of a general character, for the sake of concreteness we consider the situation where the 
wave is described by a nonlinear Schrodinger equation 

{ "ic& _e+u(r)+Mr)}0(r) = o (1) 

Here n(r) = |0(r)| 2 is the density, m is the wave mass, e is the wave energy, (3 is a constant 
and u(r) is a scattering potential, which is a random function of the coordinates. Eq.l 
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appears, for example, in the theory of electromagnetic waves propagating in nonlinear media 
[1], in the theory of hydrodynamic turbulence [2], and in the theory of turbulent plasma [3]. 

We will assume that u(r) obays the white noise statistics: (w(r)) = , (-u(r)-u(ri)) = 
l^2<5(r — ri). Here brackets () correspond to averaging over realizations of u(r), and / is the 
elastic mean free path (/ ^> k~ l = (2em)~^). In the presence of the scattering potential 
u(r), the spatial dependence of the n(r) exhibits random sample-specific fluctuations, which 
are called speckles. In the case of elastically scattering diffusive linear media this problem 
was considered in [4-6]. Below we will be interested in the statistics of n(r) in the case of 
nonlinear diffusive media. In particular, we will show that, for given boundary conditions, 
the number of solutions of Eq.l increases exponentially with sample size. A brief summary 
of some of the results has been published in [7]. 



II. NONLINEAR SPECKLES. 

Let us consider the case where a coherent wave 0o( r ) = \^no e%kr with momentum k is 
incident on a disordered sample of size L 3> / (See the insert in Fig. la). 

The r-dependence of the average density (n(r)) can be described by the diffusion equa- 
tion, 

(iM> = -D 9 -^ 

dw(j(r)} = (2) 

which is equivalent the to calculation of the diagrams shown in Fig.2.a. Here j(r) is the 
current density. In the limit A;" 1 <C / <C L the expression for the diffusion coefficient 
D = ^ has a classical form. In the case of the geometry shown in Fig.l, the total flux 
through the sample is zero, and the average density (n(r)) = no is spatially uniform. 

The term f3\(fi(r)\ 2 in Eq.l plays the role of an additional scattering potential. It will be 
shown that its contribution to the diffusion coefficient D can be neglected at small enough 
intensities of the incident beam, when 
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I« 2 « (3) 

At (3 > the propagation of a uniform wave in an infinite medium described by Eq.l is 
unstable at arbitrarily small n Q due to the effect of nonlinear self focusing. The characteristic 
length at which the self-focusing takes place is of order r( s /) ~ ek/(/3n ) 2 m [1]. Thus the 
inequality Eq.3 is equivalent to r( s /) ^> I. One can say that the regime Eq.3 corresponds to 
a system of randomly distributed weak concave and convex lenses. 

The diffusion equation approximation completely neglects interference effects, which lead 
to the existence of speckles. To describe them one should solve Eq.l before averaging over 
the realizations of -u(r). It is convenient to expand the density 

D 00 

/?n(r) = ^=E* 1/3 ^( r ) ( 4 ) 
in a complete set of ortho normal eigenstates rij(r) of the diffusion equation (/ drn 2 {v) = 1): 

where Ei ~ j?i 2 ^ are the eigenvalues of Eq.5, and i = 1,2... labels the eigenstates. We 
will show below why the expansion Eq.4 is convenient. We assume boundary conditions 
corresponding to zero current through a closed boundary, and n(r) = at the open boundary. 

Let us first substitute Eq.4 into the nonlinear term of Eq.l, and regard the = 
Ui, ...Uj, ... as independent parameters. This gives a linear equation for 0(r) 

/ 1 r) 2 D 00 \ 

[-2^9? ~ e + M W + E ^Mr) j <Kr) = (6) 

Denoting a solution of Eq.6 at a given set of parameters {ui}, as 0(r, (-u(r)}, {-u m }) and con- 
structing n(r , {u(r}) , {ui}) = |0(r, {u(r)}, {ui})\ 2 , we can write the self-consistency equa- 
tions for {Ui}: 

= m*})- (7) 

Here 
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TH^KjP (8) 

and 

= J dm(v, {u(v)}, MK(r) (9) 

are dimensionless random functions of {-Uj}, the form of which depends on u(r). 

To investigate the properties of the solutions of Eq.7, we have to know the statistical 
properties of the random functions Fi({ui}). We can infer these properties from values of 
different correlation functions of Fi({ui}) obtained by averaging over realizations of u(r). It 
is important that the statistical analysis of the random functions Fi({ui}) is equivalent to 
the analysis of linear speckles, which has been done in [4-6]. To characterize Fi({iii}) we 
calculate the following correlation functions 

(^({^})^({^})) = S ij: (10) 
([F t ({u t + AUi}) - F({^})] 2 > ~ (Au n f at Au, ~ 1. (11) 

( g x g ) . ((r/8) v3 + (8/r) v3 ) __L_ (12) 

Here 5Fi = Fi — (Fj) and (Fi({ui})) = const which is independent of {Ui}. We present the 
derivation of Eqs. 10-12 in the next section. The simple form of Eqs. 10-12 is a consequency 
of the choice of rij(r) in Eqs. 4, 5. Eq.10 indicates that mesoscopic fluctuations of different 
functions F t are uncorrelated. According to Eq.12, the correlation of derivatives of Fi({ui}) 
over different u r is small for \r — s\ > 1. It will be shown in the next section that these facts 
are consequences of the choice of rij(r) in Eq.4 as eigenfunctions of Eq.5. The introduction 
of the coefficients ^^D/y/X in Eq.4 ensures Eqs.10,11 appears as above. 

Thus we arrive at the following picture: the fluctuate randomly as functions of 

{Ui} near their average, which is independent of {Hi}. According to Eq.ll, the characteristic 
period of the fluctuations is of order one. The fluctuations both of different functions, and 
the same functions with respect to different Ui are uncorrelated. 



Using this information about the Fi({ui}) we can estimate the number of solutions of 
Eq.7 (or Eq.l). If 7 <C 1, Eq.7 has a unique solution while for 7 ^> 1 Eq.7 has many 
solutions. Let us consider the i th equation in Eq.7 and fix all variables Uj other than Mj. 
Then, at i 2 ^^ 1 3> 1 the equation has a unique solution, while at i 2 / 3 7~ 1 <C 1 the number 
of the solutions is of order r )i~ 2 l z . In Fig. la we show a qualitative graphical solution of Eq.7, 
which corresponds to the intersection of two functions: Ui, ....) and 7 _1 i 2//3 -Uj. The 

solutions are distributed in an interval of order r )i~ 2 l z . 

Therefore, to estimate number of solutions of Eqs.7 at 7 > 1, we have to take into 
account only a subset of Eqs.7 with i < I = 7 s / 2 . Since both the amplitude of fluctuations 
and the periods in the i th direction of randomly rippled hypersurfaces Fj({-Uj}) are of order 
unity, the number of solutions AT of Eqs.1,7 is proportional to the volume of an /-dimensional 
hyperparallelepiped with sides of order 7^ 2 / 3 , i < I. As a result we have 

A r ~7 / II^ 2/3 ~exp(-7 3 / 2 ) (13) 
1 ^ 

Thus the number of the solutions N of Eq.l increases exponentially with the sample size L. 

To illustrate Eq.13 we consider the case / = 2 (7 ~ 1). Then Eqs.7 can be viewed 
as two surfaces z = Fi(ui, U2) and z = F 2 (ui, U2), which are intersected by two planes 
z = 7~ 1 Mi and z = 7~ 1 2 2 / 3 w 2 respectively. A result of these intersections is two systems of 
lines in the plane u\, U2 shown in Fig. lb. The solid lines correspond to intersections between 
z = Fi(ui,u 2 ) and z = 7 -1 Wi, and are located within a strip in the U\ direction of width 
of order 7. The dashed lines correspond to the intersection of the surfaces z = F 2 (ui,u 2 ) 
and z = 7 _1 2 2 / 3 M2, and are located within a strip of width 72~ 2 / 3 in u 2 direction. The 
intersections of solid and dashed lines in Fig. lb correspond to solutions of the system of 
Eqs.7. According to Eqs. 10-12, the dashed and the solid lines in Fig. lb are uncorrelated. 
The typical distance between, say, solid lines is of order one. As a result, the number of 
solutions Nj=2 of Eqs.7 in this case is of the order of the area of a parallelogram with sides 
7 and 2" 2 / 3 7 

N I=2 ~ 7 x 2" 2/3 7. (14) 
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Now let us estimate the corrections to the diffusion coefficients originating from scattering 
from the potential (3n(r). In the case |r — r'| < I we have (D = 3) [6] 

W r W r ')> = ( 15 ) 

independently of the values of the parameters {Hi}. In the Born approximation, the nonlinear 
mean free path corresponding to scattering from the potential (5n{r) is: 

/( - ;) = (16) 

Thus the criterion Eq.3 is equivalent to Z( n ;) = r( s /) I. 

In the opposite limit /(„;) <C /, the scattering mean free path is determined by the 
scattering from the potential /3n(r), and one can neglect the random potential u(r) in Eq.l. 
In this case one can estimate the number of solutions of Eq.l by substituting instead of 
/ into Eqs.8, or by substituting 7^) into Eq.13 instead of 7, where 



3 (n (3f ( Lm ^ /2 ^ 



1{nl) 2 e v ek 

In the conclusion of this section we would like to discuss the condition 7 > 1 for the exis- 
tence of multiple solutions of Eq.l. In the absence of the nonlinear term (3n(r) the solution 
of Eq.l corresponds to particles traveling along diffusion trajectories. In the presence of the 
nonlinear term f3n(r), the probability amplitude for traveling along a diffusive trajectory ac- 
quires an additional phase of order 5x( n i) — 2 f ds5n(r) ~ 1, where the integration is taken 
along a diffusive trajectory. To estimate the value of the additional phase let us calculate 
the integral 

((>X(no) 2 ) = [jf I dsds* (5n{v)5n{v>)) (18) 

The correlation function of densities is given by Eq.15 and the estimate for Eq.18 is ((3k) 2 ^-. 
Thus the criterion 7 > 1 corresponds to {(5x(ni)) 2 ) > 1- Another interpretation is that at 
7 > 1 the sensitivity of the solutions of Eq.l to a change of u(r) increases significantly as 
compared to its single particle value [7]. 
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III. CALCULATION OF CORRELATION FUNCTIONS OF EQS.11-13 



Now we turn to the calculation of correlation functions Eqs. 10-12. According to the 
definition Eq.9 we get 

(Fi({ui})Fj({ui + Aui})) = J 2 / drdr'n i (r)n j (r')(n(r,{u l })n(r 1 ,{u l + Au,})> (19) 

Li TLq J 

In the approximation Eq.3, the value of (n(r, {■u i })n(r 1 , {u\ + Aiii})) is independent of {«/} 
and has the same form as in the linear case [4,6]. Therefore, one can use the standard 
diagram technique for averaging over u(r) [11]. The diagrams describing the correlation 
functions (n(r, {ui})n(r', {ui + Au t })) are shown in Fig. 2b, c. Alternatively, one can solve 
the Langevin equation [4,6], valid at |r — r'| ^> I 

divSj(r) = (20) 

5j(r) = - B S Mr) + jL(r ' Wr)} ' { ^ }) - (21) 
The correlation function of random Langevin forces J L (r) is given by the diagram shown in 
Fig. 2b. It can be written as 

(J*(r, Mr}, {«,}) J/(r', {u(r)}, + A*})) = 
27r/ 

^|(0(r, {«,})^(r, {«, + A^}))| 2 ^5(r - r') (22) 

The expression for (0(r, {-u;})0*(r, {-u; + A-u;})) is given by the ladder diagrams shown in 
Fig. 2b, where the inner Green function correspond to {ui}, while the outer Green functions 
correspond to {u m + Au m }. Eq.22 is a generalization of corresponding equations introduced 
in [4,6]. Namely, in the case Aui = one has to substitute (n) 2 = n 2 , in Eq.22 instead of 
|(0(r,{^})0*(r, {u l + Au l }))\ 2 . 
Using Eqs.20-22 we get 

(n{r,{u l })n{r',{u l + Aui})) = 

Mj^ dn^)^) ^ {S(}¥ , (r _ {fii + A . }))|2 (23) 



n (r y HS ^M (24) 

I E l 

Here II(r, r') is the Green's function of Eq.5. 

At Aui = 0, doing the integral over ri in Eq.23 by parts, using Eq.22 and taking into 
account the orthogonality of the functions nj(r), we get Eq.10. The diagrams in Fig.2.b for 
(0(r, {ui})<ft*(r, {ui + Aui})) are equivalent to the equation [9,10] 

( jD £ + ^^' V3A ^ (r) ) W r '{^})^( r '{^ + A ^})) = ° (25) 

The existence of the term proportional to Aui in Eq.25 reflects the fact that the inner lines in 
the diagrams shown in Fig. 2b correspond to the Green's functions in the scattering potential 
characterized by while the outer lines correspond to those of {v,i + Aw;}. Note that as 
long as Eq.3 holds Eq.25 is independent of 

Solving Eq.25 by perturbation theory with respect to Aui we get 

(0(r, {u l })<j>*(r,{u l + Au l })) = 

no ^ +i ^ t ^ AUt _ el / w? *»(y>>% ...) . (26) 

We can neglect the second term in brackets in Eq.26 because it is of order (i~ 1 / 3 A'Uj) (and 
its contribution to Eq.22 is of order i~ 2 ^ 3 (Aui) 2 ), while the contribution to Eq.22 from the 
third term is (Aiti) 2 . To get the latter estimate we took into account that n(r, r') ~ (DL)^ 1 
at |r - r'| ~ L. Substituting Eq.26 into Eq.19,22,23, we get Eqs.11,12. 

IV. DISCUSSION 

The estimate Eq.14 was made for the case of a typical realization of the scattering 
potential. On the other hand, even at 7 < 1, there are rare realizations of u(r), which 
correspond to several solutions of Eq.l. 

The results presented above hold for arbitrary sign of (3. This is quite different from the 
situation in the pure case (u(r) = 0) [1] where at (3 > self-focusing takes place. 
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At 7 3> 1, the solutions of Eq.l exhibit exponentially large sensitivity to changes of 
parameters of the system [7]. Consider, for example, the case where the incident angle 9 of 
the wave is changing, and suppose that a solution of Eq.l follows this change adiabatically. 
Then an exponentially small change 

A6 ~ exp(-jj 7 3 / 2 ) (27) 

will lead to disappearance of the solution, and the system will jump to another solution. 

Similar phenomenon may occur in the system of interacting electrons is disordered metals: 
it can be unstable with respect to the creation of random magnetic moments. This would 
correspond to Finkelshtain's scenario [12]. However, in this case in order to get a self- 
consistency equation, which would be an analog of Eq.7, we have to integrate over electron 
energies, which decreases the amplitude of mesoscopic fluctuations. As a result, the situation 
with many solutions may occur only in the D=2 case and the characteristic spatial scale will 
be of the order of the electron localization length in the linear problem. Thus the problem 
of interacting electrons in disordered metals remains unsolved. 

We would like to mention a similarity of the problem considered above to the problem 
of spin glasses. To illustrate this point let us consider a model in which the coefficient /5(r) 
in Eq.l is nonzero only at points r = r a , a — 1, 2 

/3(r)=/3 ££/(r-r a ), (28) 

a 

where t/(|r|) is a short range function decaying on characteristic distance R < 1/k, and of 
a maximum height Uo , and are randomly distributed in space with given density. Then 
Eq.l can be rewritten only in terms of values of (fi a = 0(r = r a ) 

<p a = foUoR 3 E G(r a , v p ) (29) 

P 



G(r a , Tp) ~ " v ' | _ m 1/2 v P> |r a - r fi \ » I (30) 

where G(r a , rp) is the Green function of the linear Schrodinger equation, and the phase 
S(r a ,rp) is a random quantity at |r a — rp\ ^> I. The major difference between Eq.30 and 
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the spin glass problem is that in the former case one is interested in the minimum of the 
free energy, while in the case of Eq.30 the boundary conditions are given. Thus there are no 
thermodynamic criteria on how to choose between multiple solutions of the stationary Eq.l. 
One of the possibilities is that the state of the system is determined by its history. 

We may question how many of these stationary states are stable. The time-dependent 
non-linear Schrodinger equation may be obtained from Eq.l by substituting e for id t . Equiv- 
alently, we can write equations for Ui(t). We still assume the same stationary boundary 
conditions. In the absence of a complete solution we present a qualitative picture. 

The characteristic time of change of the i'th harmonic Ui is r% ~ E~ x . Thus we have 

t % (1 + 7&({«i})) dMt) = 7 rlF({^}) - Hi (31) 

Here #«({««}) ~ 1 is a dimensionless random function. The statistical properties of the 
function ^ ({««}) are roughly the same as those for Fi({ui}). Namely, the amplitude and 
sign randomly oscillate with a characteristic periods of order one. Strictly speaking, Eq.31 
holds if \rid t Ui\ <C 1 and the characteristic time of establishing stationary distributions n(r) 
at given Ui(t) is much shorter than the characteristic time of change of Ui{t). In other words, 
Eqs.31 represent a sort of hydrodynamics. This takes place, for example, near the critical 
instability points (see, for example, the point "a" on Figl.a.). Generally, Tid t Ui(t) ~ 1, so 
we have to keep not only the first but also several higher time derivatives in Eqs.31. We 
believe, however, that the model Eqs.31 captures the qualitative features of the system's 
dynamics correctly even in this case. 

Linearizing Eqs.31 near the stationary solutions, we arrive at the conclusion that the 
fraction of solutions of Eq.7 which are stable is of order 2 _/ . Thus, at 7 3> 1 the number of 
stable stationary solutions is still exponentially large. 

In principle, Eqs.31 can also have nonstationary solutions as t — > 00. Obviously, the 
characteristic amplitudes of the solutions are given by Eq.7, (\5ui(t)\ < 7?~ 2 / 3 ). A complete 
investigation of {-Uj(t)} is a complicated and still unsolved problem. For example, for 1=1 
and t — > 00 only stable stationary solutions are relevant. For 1 = 2, and t — > 00, depending 
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on the properties of the realizations of the random sample specific functions Fi(ui,u 2 ) and 
F 2 (-Ui,m 2 ), one can, additionally, have periodic, quasiperiodic, and chaoric in time solutions. 
If I = 3, one can also, have strange (or stochastic) attractors as solutions of three differential 
equations. 

In this respect we would like to mention papers [13-15], where an attempt to describe 
temporal nonlinear speckles was done. We believe that these results are incorrect. To 
estimate the instability threshold it is sufficient to expand the nonlinear Langevin equations 
in powers of (3 (see diagrams shown in Fig. 2 in [7], or in Fig. 4 in [14]). This is why the 
authors of [13-15] were able to reproduce the instability criterion 7 > 1 obtained in [7]. 
Strictly speaking, this approach holds only in the case when 7 <C 1. For 7 > 1 this approach 
describes the system incorrectly. In the absence of solid mathematical procedure the authors 
of [13-15] made assumptions about the nature of solutions of Eq.l beyond the instability 
point 7 > 1. Namely, they made assumptions that at t — > 00 the finction n(r,t) must 
exhibit oscillations in time and they assumed some form of time correlations of n(r, t). Both 
assumptions are incorrect. This can be seen, for example, from the fact that their approach 
cannot reproduce the existence of an exponentially big (for 7 ^> 1) number of multiple 
stationary solutions of Eqs.7, which are singular points of Eqs.31. 

The simplest situation where the deficiency of the approach of [13-15] is most evident 
takes place near the first instability point 7 = 7 C ~ 1. This is schematically shown by the 
dashed line in Fig. la. The critical value 7 C is sample specific. In this case, say, the first 
equation (i — 1) in the system of Eqs.7 has three solutions. Two of them are close to each 
other (see the point "a" in Fig. la.). Let us start with a discussion of the system dynamics 
near this point. We would like to stress that in this case {T\d t ui{t)) / u\ <C 1 and our analysis 
is rigorous. This is exactly the regime considered in [15]. Since d t Ui/ui ~ E { (7 — 7 c )£'i 
for i > 1, we can neglect the time derivatives in all equations of the system Eqs.31, except 
the one with % = 1. Moreover, since this point of instability is a rare one, all equations 
in Eqs.7 with i > 1 have unique solutions. Thus, the system's dynamics is described by 
just one first order differential equation for u\. Near the point "a" in Fig.la the functions 
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gi({ui}) and Fj>i({-Uj}) change slowly. Thus, one stationary solution is stable and the other is 
unstable. The solution indicated by the letter "b" in Fig. la is also stable because it is related 
adiabatically to the unique stable solution for 7 < j c . At t — > 00 the system approaches one 
of the stable stationary solutions, described by Eqs.7 [7]. Thus the assumption made in [15] 
about existence of oscillating in time solutions near the critical point is wrong. 

Far from the critical point, when I ^> 1 and the number of relevant equations of the 
system of Eqs.31 is larger than one, depending on initial conditions and the form of random 
functions Fj({iZj}, Eqs.31 can have periodic, chaotic solutions and strange attractors, in 
addition to stationary points. The fractions of the phase space which at t — > 00 are attracted 
to these types of motion are currently not known. In any case, the solutions of the time- 
dependent nonlinear Schrodinger equation, or Eqs.31, have a very complicated non-Gaussian 
character, which is very different from that assumed in [13-15]. 

Finally, we would like to mention that in reality, in the case of very large N, the time 
dependent fluctuations of external sources become important. They are not described by 
Eq.l and the problem remains unsolved. 
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FIGURES 
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FIG. 1. a) Graphical solution of Eq.7. The wavy line corresponds to and the straight 

lines corresponds to 7~ 1 i 2 / 3 Uj at different values of 7. The dashed line illustrates the critical 
instability point, when a solution of Eq.7 becomes nonunique. b) The solid lines correspond to 
the intersection F\(u\,U2) and 7~ 1 ui, while the dashed lines correspond to the intersection of 
F 2 (ui,v,2) and 7~ 1 2 2 / 3 u 2 . 
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FIG. 2. a) Diagrams describing (n(r)). Solid lines correspond to Green functions of Eq.l with 
(3 = 0. Dashed lines correspond to 7r<5(r — r')/lm 2 . b) and c) Diagrams describing Eq.24. The 
inner solid lines describe the Green functions which correspond to {ili}, while the outer solid lines 
correspond to {u« + Auj}. 
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